PERFECT SIMULATION OF VERVAAT PERPETUITIES 



JAMES ALLEN FILL AND MARK HUBER 

Abstract. We use coupling into and from the past to sample perfectly 
in a simple and provably fast fashion from the Vervaat family of perpetu- 
ities. The family includes the Dickman distribution, which arises both 
in number theory and in the analysis of the Quickselect algorithm, 
which was the motivation for our work. 



1. Introduction, background, and motivation 

1.1. Perpetuities in general. Define a perpetuity to be a random vari- 
able Y such that 

Y = Wi + WiW2 + WiW2W3 + --- (1.1) 

for some sequence Wi, W2, W3, ... of independent and identically distributed 
random variables distributed as W. Throughout this paper we assume 
W >0 (a.s.) and El^ < 1 since these simplifying assumptions are met 
by the Vervaat perpetuities, which are discussed in Section [1.21 and are the 
focus of our work. [Some authors define a perpetuity as in (jl.ip but with 1 
added to the right-hand side.] The distribution of such a random variable Y 
is also referred to as a perpetuity. General background on perpetuities is 
provided in the first paragraph of Devroye [U Section 1]. To avoid repe- 
tition, we refer the reader to that paper, which also cites literature about 
perpetuities in general and about approximate sampling algorithms. 

Within the general framework of (jl.ip . the following simple observations 
c 

(with = denoting equality in law, or distribution) can be made: 

(i) The random variable y > is finite almost surely; indeed, its ex- 
pected value is finite: 

BY = < 00. 

l-EW 
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(ii) The perpetuity satisfies the distributional fixed-point equation 

Y = W{1 + Y) (1.2) 

where, on the right, W and Y are independent. [In fact, this fixed- 
point equation characterizes the distribution (11. ip .] 

(iii) If W has a density fw, then Y has a density /y satisfying the 
integral equation 

fY{y)= ri^ + vr'fw(-^)fYiv)dv, y>o. 

Jo Kl + r]/ 

1.2. Quickselect, the Dickman distribution, and the Vervaat fam- 
ily of perpetuities. Our interest in perpetuities originated with study 
of the running time of the Quickselect algorithm (also known as Find), 
due to Hoare [1]. Quickselect(n, m) is a recursive algorithm to find the 
item of rank m > 1 (say, from the bottom) in a list of n > m distinct 
numbers (called "keys"). First, a "pivot" is chosen uniformly at random 
from among the n keys and every key is compared to it, thereby deter- 
mining its rank (say, j) and separating the other keys into two groups. 
If j = m, then Quickselect(n, m) returns the pivot. If j > m, then 
Quickselect(j — l,m) is applied to the keys smaller than the pivot. If 
j < m, then Quickselect(n — j^m — j) is applied to the keys larger than 
the pivot. 

Let C(n, m) denote the (random) number of key comparisons required by 
the call Quickselect (n, m), and write l(^) to mean the indicator that has 
value 1 if the Boolean expression A is true and otherwise. Let Jn denote 
the rank of the first pivot chosen. Immediately from the description of the 
algorithm we find the distributional recurrence relation 

C(n, m) = (n — 1) + 1( J„, > m)C( J„ — 1, m) + 1( Jn < m)C*{n — Jn, m — Jn), 

(1.3) 

where, on the right, (i) for each fixed r and s the random variable C(r, s) is 
distributed as the number of key comparisons required by Quickselect (r, s), 
the joint distribution of such random variables being irrelevant; (ii) similarly 
for C*{r, s); (iii) the collection of random variables C(r, s) is independent of 
the collection of C*(r, s); and (iv) J„ is uniformly distributed on {1, . . . , n} 
and is independent of all the C's and C*'s. 

The distribution of C(n, m) is not known in closed form for general finite n 
and m, so we turn to asymptotics. For any fixed m, formal passage to the 
limit as n — )■ oo suggests that 

C{n,m) 

Y [n, m) := 1 

n 

has a limiting distribution C{Y) satisfying the fixed-point equation 



Y = U{1 + Y). 



PERFECT SIMULATION OF VERVAAT PERPETUITIES 



3 



This is indeed the case, as was shown by Mahmoud et al. [12]. Recalhng 
the characterization (]1.2p . we see that the law of 1" is a perpetuity, known 
as the Dickman distribution. Curiously, if m is chosen uniformly at random 
from {1, . . . , n} before the algorithm is run, then the limiting distribution of 
Y(n,m) is the convolution square of the Dickman distribution, as was also 
shown in |12j . 

See [HI Section 2] concerning various important settings, including number 
theory (largest prime factors) and combinatorics (longest cycles in permuta- 
tions), in which the Dickman distribution arises; also see [6] for some basic 
facts about this distribution. The support of the distribution is [0, oo), and 
simple integral expressions are known for the characteristic and moment- 
generating functions of this log-concave (that is, strongly unimodal) and 
infinitely divisible distribution. In addition, the Dickman distribution has 
a continuous density / that is constant with value e~'^ over (0, 1] (here 7 is 
Euler's constant). Over each interval {k,k + 1] with k a positive integer, / 
is the unique solution to the delayed differential equation 

f'iy) = -fiy - i)/y. 

Still, no closed form for / is known. 

Some years back, Luc Devroye challenged the first author to find a method 
for simulating perfectly from the Dickman distribution in finite time, where 
one assumes that only perfect draws from the uniform distribution and ba- 
sic arithmetic operations such as multiplication (with perfect precision) are 
possible; in particular, numerical integration is not allowed. The problem 
was solved even prior to the writing of Devroye's paper [T], but the second 
author pointed out extensive simplification that could be achieved. The re- 
sult is the present collaborative effort, which shows how to sample perfectly 
in a simple, provably fast fashion from the Dickman distribution, and more 
generally from any member of the Vervaat [15] family of perpetuities han- 
dled by Devroye [1]. To obtain the Vervaat perpetuities, one for each value 
of < /3 < 00, choose W = C/^/^ in (fr2D . 

1.3. A quick review of Devroye's method. To compare with the Markov- 
chain-based method employed in the present paper, we first review high- 
lights of Devroye's [1] method for perfect simulation of perpetuities. Devroye 
sketches a general approach and carries it out successfully for the Vervaat 
family. 

The underlying idea of Devroye's approach is simple acceptance-rejection: 
1. Find explicit h > and 1 < c < 00 so that 



fy 1^ h everywhere and 




2. Generate X with density c h. 

3. Accept X with probability fy{X) /h{X)] otherwise, reject X and 
independently repeat steps 2-3. 
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Here is how Devroye [T] carries out these three steps in the case of Vervaat 
perpetuities. (It is somewhat surprising how simple the first two steps turn 
out to be.) 

1. One can take 

h{x) ■.= mm\/3{l3 + l)x-^,Pxl^-^] and c = (/3 + . 

2. Let Ui and U2 be independent uniform(0, 1) random variables and 
set/ 

X:= (/3 + l)i/(/3+i)[/i//'[/-i. 

3. Since /y can't be computed exactly, having generated X = x and 
an independent uniform random variable f/3, one must figure out 
how to determine whether or not U3 < fY{x) /h{x). 

Devroye's solution for step 3 is to find explicitly computable approxima- 
tions fn{x) to /y(x) and explicitly computable bounds Rn{x) so that, for 
every x, 

\fn{x) - fY{x)\ < Rn{x), Rn{x) ^ as rn> CX). 
Devroye's functions /„ are quite complicated, involving (among other things) 

• the sine integral function S{t) := Jq^^j^ ds and approximations to 
it computable in finite time; 

• explicit computation of the characteristic function (pY of Y; 

• use of quadrature (trapezoidal rule, Simpson's rule) to approximate 
the density fy as the inverse Fourier transform of 0y. 

Devroye proves that the running time of his algorithm is finite almost 
surely (for any /3 > 0), but cannot get finite expected running time for any 
/3 > without somewhat sophisticated improvements to his basic algorithm. 
He ultimately develops an algorithm so that 

ET < oo for /3 > 5/3. 

More careful analysis would be difficult. Devroye makes no claim "that these 
methods are inherently practical". We do not mean to criticize; after all, 
as Devroye points out, his paper [1] is useful in demonstrating that perfect 
simulation of perpetuities (in finite time) is possible. 

The approach we will take is very simple conceptually, very easy to code, 
and (at least for Vervaat pereptuities with /3 not too large) very fast (prov- 
ably so). While we have hopes that our methodology can be generalized to 
apply to any perpetuity, in this paper we develop the details for Vervaat 
perpetuities for any /3 > 0. 

1.4. Our approach. Briefly, our approach to perfect simulation of a per- 
petuity from the Vervaat family is to use the Markov-chain-based perfect 
sampling algorithm known as coupling into and from the past (CIAFTP) 
([HI [9]) that produces draws exactly from the stationary distribution of a 
Markov chain. CIAFTP requires use of a so-called dominating chain, which 
for us will be simple random walk with negative drift on a set of the form 
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{xq — l,xo,xo + 1,xq + 2, . . .}, where xq is a fixed real number at least 2. In 
order to handle the continuous state space, multigamma coupling |13j will 
be employed. 

As we shall see, all of the Markov chain steps involved are very easy 
to simulate, and the expected number of steps can be explicitly bounded. 
For example, our bound is the modest constant 15 in the case /3 = 1 of 
the Dickman distribution, for which the actual expected number of steps 
appears to be a little larger than 6 (consult the end of Section [5]). 

The basic idea is simple. From the discussion in Section 11.11 it is clear 
that the kernel K given by 

K{x,-) := CiWil + x)) (1.4) 

provides a Markov chain which, for any initial distribution, converges in law 
to the desired stationary distribution C{Y), the perpetuity. 

An update function or transition rule for a Markov chain on state space 
with kernel K is a function (p : Q x 0,' (for some space O') so that for a 
random variable W with a specified distribution on we have £.{(j){x, W)) = 
K{x,-) for every x G O. There are of course many different choices of update 
function for any particular Markov chain. To employ coupling from the past 
(CFTP) for a Markov chain, an update function that is monotone suffices. 
Here, an update function (p is said to be monotone if whenever x ^ y with 
respect to a give partial order on Q we have (f){x, w) ^ <j){y, w) for all w G Q,'. 

Consider the state space [0,oo) linearly ordered by <, and let W be 
distributed as in ()1.2p . Then 

4>na.tuTal{x,w) = w{l + x) 

provides a natural monotone update function. If we wish to do perfect 
simulation, it would appear at first that we are ideally situated to employ 
coupling from the past (CFTP). 

However, there are two major problems: 

(i) The rule i;^'naturai is strictly monotone, and thus no two trajectories 
begun at distinct states will ever coalesce. 

(ii) The state space has no top element. 

It is well known how to overcome these difficulties: 

(i) Instead of ^natural) we use a multigamma coupler [13] . 

(ii) Instead of CFTP, we use CIAFTP with a dominating chain which 
provides a sort of "stochastic process top" to the state space [HJ [9] . 

Our work presents a multigamma coupler for perpetuities that is monotone. 
Monotonicity greatly simplifies the use of CFTP [T3] and CIAFTP. Useful 
monotone couplers can be difficult to find for interesting problems on con- 
tinuous state spaces; two very different examples of the successful use of 
monotone couplers on such spaces can be found in [16] and [5]. 

1.5. Outline. In Section [2] we describe our multigamma coupler; in Sec- 
tion [31 our dominating chain. Section |4] puts everything together and gives 
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a complete description of the algorithm, and Section [5] is devoted to bound- 
ing the running time. Section [6] briefly discusses approaches similar to ours 
carried out by two other pairs of authors. 



2. The multigamma coupler 

The multigamma coupler of Murdoch and Green |13j is an extension of 
the 7 coupler described in Lindvall [11] that couples a single pair of random 
variables. An update function can be thought of as coupling an uncountable 
number of random variables simultaneously, thus the need for "multi" gamma 
coupling. The goal of multigamma coupling is to create an update function 
whose range is but a single element with positive probability, in order to use 
CIAFTP as described in Section H 

Originally multigamma coupling was described in terms of densities; here 
we consider the method from the point of view of the update function. 
Suppose that the update function can be written in the form 

(j){x,w) = l{w e A;j;)(j)i{w) + l{w ^ Ax)(j)2ix,w). 

Now let A := Ox Ax- liW G A, then the Markov chain will transition from x 
to (p{x, W) = (j)i{W), for all values of x. 

Recall our Markov kernel ()1.4p . In the Vervaat case W = U^/l^ the con- 
ditional distribution of W{1 + x) given W < 1/(1 -|- x) is the unconditional 
distribution of W. 

Proposition 2.1. For any Vervaat perpetuity, the function 

<l)ix,wil),w{2)) := 1 (^wil) < ^) w{2) + l (^w{l) > w{l){l + x) 

is a monotone update function for the kernel K at (|1.4p . More explicitly, 
(j){x,w{l),w(2)) is nondecreasing in x G [0,oo) for each fixed w{l),w{2) € 
[0,1), and the distributions of (l){x,W{l),W{2)) andW{l + x) are the same 
when W{1), W{2) are two independent random variables distributed asW = 
U^^^, where U is uniformly distributed on [0,1). 

Proof. Since the conditional distribution of U{1 + x)^ given U < [1 + x)~^ 
is uniform on [0, 1), by taking /3th roots we find 

C{W{1 + x)\W <{l + x)-^) = C{W), 

independent of x, as remarked prior to the statement of the proposition. 
Therefore (j){x, W{1), W{2)) has the same distribution as W{\ -\- x). 

Also, for fixed w{l),w{2) € [0,1) the function (f){-,w{l),w{2)) has the 
constant value w{2) < 1 over [0,w{l)~^ — 1] and increases linearly over 
[w{l)~^ — l,oo) with value 1 at the left endpoint and slope u'(l), so mono- 
tonicity is also clear. □ 
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3. A DOMINATING CHAIN 

Dominating chains were introduced by Kendall [8] (and extended in [9]) 
to extend the use of CFTP to chains on a partially ordered state space with 
a bottom element but no top element. A dominating chain for a Markov 
chain (Xt) is another Markov chain (Dt) that can be coupled with (Xt) so 
that Xt < Dt for all t. In this section we give such a dominating chain for the 
Vervaat perpetuity, and in the next section we describe how the dominated 
chain can be used with CIAFTP. 

Since our multigamma coupler (j) of Proposition 12.11 is monotone, if we 
can find an update function ijj{x,w{l),w{2)) on a subset S of the positive 
integers x such that 

(j){x,w{l),w{2)) <'ip{x,w{l),w{2)) for all xG 5, (3.1) 

then [with the same driving variables W{1),W{2) as in Proposition 12. Ij ^ 
drives a dominating chain. This is because if x G [0, oo) and y G S" and 
X < y, then 

(Pix,w{l),wi2)) < <P{y,wil),wi2)) < ij{y,w{l),w{2)). 

We exhibit such a rule in our next result, Proposition 13. 1[ It is im- 
mediate from the definition of ■0 that the dominating chain is just a simple 
random walk on the integers {xq — 1, xq, . . . } that moves left with probability 
2/3 and right with probability 1/3; the walk holds at — 1 when a move to 
xq — 2 is proposed. In the definition of ■0, note that no use is made of w(2). 

Proposition 3.1. Fix < (3 < oo. Let (p be the multigamma coupler for 
Vervaat perpetuities described in Proposition \2.1\ Define 

2 



Xq 



1 > 2 (3.2) 



1 - (2/3)i//3 
and, for x ^ S := {xq — 1, xq, . . .}, 

'ilj{x,w{l),w{2)) ■.= x + l(w(l) > (2/3)^/^) - l(w(l) < (2/3)^/^, x > xq). 

Then (13. ip holds, and so ip drives a chain that dominates the (p-chain. 

Proof. To establish ()3.ip . suppose first that x > xq. Since 

(j){x,w{l),w{2)) <x + l and ^{x,w{l),w{2)) e {x - l,x + 1}, 

we need only show that if ^p{x,w{l),w{2)) = x — 1 [i.e., if w{l) < (2/3)^/^^], 
then (p{x, w{l),w{2)) < x — 1. Indeed, if w{l) < 1/(1 + x), then 

(/)(x, w{l),w{2))) = w{2) < 1 < xo - 1 < X - 1; 

and if w{l) > 1/(1 + x), then 

(l){x,w{l),w{2)) ='u;(1)(1 + x) < (2/3)^/^ (x + 1) < x - 1, 

where the last inequality holds because 

x-1 2 2 /2y^^ , 

= 1 > 1 > - • (3.3 

x + 1 x + l~ xo + l~V3y 
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[This is how the value of xq was chosen in ()3.2p .] 

The proof for x = xq — 1 is quite similar, so we omit most details. In 
place of (I3.3P one needs to establish 



but this follows from the last inequality in (13. Sp since xq > 2 > 1. 

4. A PERFECT SAMPLING ALGORITHM FOR VeRVAAT PERPETUITIES 

For an update function (^(x, w) and random variables W-t, . . . , VF-i, set 



so that if X-t = X, then Xq = F^^{x). To use coupling from the past with 
a dominated chain for perfect simulation, we require: 

(i) A bottom element for the partially ordered state space of the 
underlying chain X, a dominating chain D, and update functions 
(j){x, w) and ip{x^ w) for simulating the underlying chain and domi- 
nating chain forward in time. 

(ii) The ability to generate a variate from the stationary distribution of 
the dominating chain. This variate is used as Dq, the state of the 
dominating chain at time 0. 

(iii) The ability to simulate the dominating chain backwards in time, so 
that -D-i, -0-2, • • • , D^t can be generated given Dq. 

(iv) The ability to generate i.i.d. random variables W-t, W-t+\, . . . , W-i 
for the update function conditioned on the values of -Dqj • • • j D-f 

(v) The ability (given W-t, ■ ■ ■ , W-i) to determine whether F^^{x) takes 
on a single value for all x G [0,-D_j]. (This detection of coalescence 
can be conservative, but we must never claim to detect coalescence 
when none occurs.) 

With these pieces in place, dominated coupling from the past [9] is: 

1. Generate Dq using (ii). 

2. Let t' ^l,t^O. 

3. Generate D_t-i, . . . ,D-t' using (iii). 

4. Generate W-f, . . . , W-t~i using D-t', ■ ■ ■ , D-t and (iv). 

5. If F^j/([0, D_t']) = {x}, then output state x as the random variate 
and quit. 

6. Else let t ^ t' , t' ^ t' + 1, and go to step 3. 

The update to t' in step 6 can be done in several ways. For instance, 
doubling time rather than incrementing it additively is often done. For our 
chain step 5 only requires checking to see whether W-f < 1/(1-|-L'_j'), and 
so requires constant time for failure and time t' for success. Thus step 6 as 
written is efficient. 




F%{x) = (/>((/>(. . . (t>{x, W-t), . . . , W-2), W-i), 
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Because the dominating chain starts at time and is simulated into the 
past, and then the underlying chain is (conditionally) simulated forward, 
Wilson [16] referred to this algorithm as coupling into and from the past 
(CIAFTP). 

Now we develop each of the requirements (i)~(v) for Vervaat perpetuities. 
Requirement (i) was dealt with in Section [21 where it was noted that the 
dominating chain is a simple random walk with probability 2/3 of moving 
down one unit towards and probability 1/3 of moving up one unit. The 
chain has a partially absorbing barrier at — 1, so that if the chain is at 
state xo — 1 and tries to move down it stays in the same position. 

The stationary distribution of the dominating random walk can be ex- 
plicitly calculated and is a shifted geometric random distribution; require- 
ment (ii) can be satisfied by letting Dq be xq — 2+G where G € {1, 2, ... } has 
the Geometric distribution with success probability 1/2. [Recall that we can 
generate G from a Unif(0, 1) random variable U using G = [—(In U) / (In 2)] .] 
Simulating the D-chain backwards in time [requirement (iii)] is also easy 
since D is a birth-and-death chain and so is reversible. 

Now consider requirement (iv). Suppose that a one-step transition D^t+i 
= D-t + 1 (say) forward in time is observed, and consider the condi- 
tional distribution of the driving variable W-t that would produce this [via 
D^t+i = V'(-D-t) W^-t)]- What we observe is precisely that the random 
variable VF_t(l) = U]/f fed to the update function ^p must have satisfied 
W-t{l) > (2/3)1/^, i.e., U^t > 2/3. Hence we simply generate W-t{l) 
from the distribution of u]/f conditionally given U-t > 2/3. Similarly, if 
D^t+i = D^t — 1 or D^t+i = D^t = xq — 1, then we generate W-t{l) 
conditioned on W-t{l) < (2/3)^^^, i.e., on U-t < 2/3. The random variable 
W-t{2) is always generated independently as the 1//3 power of a Unif(0, 1) 
random variable. 

Finally, requirement (v) is achieved by using the multigamma coupler 
from the last section; indeed, if ever W-t{l) < l/{D-t + 1)> then the un- 
derlying chain coalesces to a single state at time — t -|- 1, and hence also at 
time 0. 

Thus our algorithm is as follows: 
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1. Generate G ~ Geom(l/2) and let Dq ^ xq-2 + G. Set t ^ 1. 

2. Having generated generate D^t by moving up from 

with probabihty 1/3 and "down" with probabihty 2/3. Here, and 
in step 4 below, a move up adds 1 and a move "down" subtracts 1 
unless -D_(t_i) = xq — 1, in which case D^t = — 1, too. 

3. Read the transition forward from D^t to D^t+i- 

4. If the forward move is up, impute U-t ~ Unif(2/3, 1); if the forward 
move is "down", impute U-t ~ Unif (0,2/3]. In either case, set 

W.t{l) = U^Jf. 

5. If W-t{l) < l/{D-t + 1), then generate W-t{2) independently and 
set X-t+i = W-t{2); otherwise increment t by 1 and return to 
step 2. 

6. Using the imputed driving variables W-s{l) {t > s > 1), and gen- 
erating independent variables W-s(2) as necessary, grow the X- 
trajectory forward to time 0: 

Here (p is the multigamma coupler described in Proposition 12.11 

7. Output Xq. 

5. A BOUND ON THE RUNNING TIME 

Theorem 5.1. Let T be the number of steps taken backwards in time. (This 
is also the number of times steps 3, 4, o-nd 5 of the algorithm at the end of 
Section^^ are executed.) Define xq = [2/(1 — (2/3)^/^)] — 1 as in Proposi- 
tion \3.1[ Then for any < /3 < cxo we have 

^ <BT <2{xo + lf + 3 (5.1) 



X 



and hence ET = e'^''^'^+®('^) as ^ ^ oo; and 

ET^las/3^0. (5.2) 

Proof. The lower bound in (15. ip is easy: Since Dt > xq — 1, the expectation 
of the number of steps t needed to get (U-t)^^^ < l/(I?_t + 1) is at least 
Xq. We proceed to derive the upper bound. 

Consider a potential function $ = (^'t)t>o such that $i depends on 
D-t, ■ ■ ■ ,Do and W-t, ■ ■ ■ , W-i in the following way: 

Jo if t > r 

* ~ \ D^t - (xo - 1) + |(xo + 1)^ otherwise. ^^"^^ 

We will show that the process (<I>tAr + (l/3)(i AT)) is a supermartingale and 
then apply the optional sampling theorem to derive the upper bound. Let 
be the cr-algebra generated by -D_t, . . . , Dq and W-t, ■ ■ ■ , W-i, so that T 
is a stopping time with respect to the filtration {Tt). Suppose that T > t. 
In this case there are two components to ^t — ^t-i- The first component 
comes from the change from D-t to D-t+i- The second component comes 
from the possibility of coalescence (which gives T = t) due to the choice of 
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W-t{l). The expected change in $ is the sum of the expected changes from 
these two sources. 

For the change in the Z)-chain we observe 

= >xo-l) + il(L»_m =^0-1). (5.4) 

The expected change in <^ due to coalescence can be bounded above by 
considering coalescence only when D-t+i = xq — 1, in which case we have 
D_t G {xq — l,xo}. But then coalescence occurs with probability at least 
[1/(1 + xq)]^ , and if coalescence occurs, drops from {2/3){xq + l)'^ down 
to 0. Therefore the expected change in <I> from coalescence when D-t+i = 
xo - 1 is at most -(2/3)(xo + lf/{xo + 1)^ = -2/3. Combining with 
yields 

E[$t-$t_il7;_i] 

< (1 _ |)l(D_t+i > - 1) + (i - l)l{D.t+i = xo - 1) = -i. 

So whenever T > t, the potential <I> decreases by at least 1/3 on average at 
each step, and hence ($jat + (l/3)(t A T)) is a supermartingale. Since it 
is also nonnegative the optional sampling theorem can be applied (see for 
example [3], p. 271) to yield 

E$r + ^ET < E$o- 

Since = 0, it remains to bound E<&o- Recall that Dq — (xq — 2) is a 
geometric random variable with parameter 1/2 and thus has mean 2. Hence 
E$o = 1 + (2/3)(xo + I)'', giving the desired upper bound on ET. 

We now prove ()5.2p . using notation such as T{/3) to indicate explicitly the 
dependence of various quantities on the value of the parameter /3. However, 
notice that xo(/3) = 2 for all < /3 < /3o := ln(3/2)/ln3 and hence that 
the same dominating chain D may be used for all such values of /3. For 
0<P <Po define 



Conditionally given the entire D-process, since 1 — {Ds + l)"'^ is the prob- 
ability that coalescence does not occur at the sth step backwards in time, 
Yt{l3) equals the probability that coalescence does not occur on any of the 
first t steps backwards in time. Thus 

P{Ti(3)>t) = BYti(3), t = l,2,..., 

and hence 

oo 

ET{(3) = l + Y,^yt{P)- (5.6) 

t=i 
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We now need only apply the dominated convergence theorem, observing 
that Yt{(3) > is increasing in < /3 < /3o and that 1 + Et^i E i^t(/3o) = 
Er(/3o) < 2 -3* + 3 = 6 < oo by dsn)- □ 

For any fixed value of /?, it is possible to find ET to any desired precision. 
Let us sketch how this is done for the Dickman distribution, where f3 = 1. 
First, we note for simplicity that T in Theorem 15 . 1 1 has the same distribution 
as the time it takes for the dominating random walk D = {Dt)t>o, run for- 
ward from the stationary distribution at time 0, to stop, where "stopping" 
is determined as follows. Having generated the walk D through time t and 
not stopped yet, let Ut be an independent uniform(0, 1) random variable. If 
Ut > 2/3, let D move up (and stopping is not possible); ii Ut < 2/3, let D 
move "down" and stop (at time t + 1) if C/^ < l/{Dt + 1). Adjoin to the 
dominating walk's state space {xq — 1, xq, xq + 1, • • •} an absorbing state cor- 
responding to stopping. Then the question becomes: What is the expected 
time to absorption? We need only compute the expected time to absorp- 
tion from each possible deterministic initial state and then average with 
respect to the shifted-geometric stationary distribution for D. For finite- 
state chains, such calculations are straightforward (see for example [7], pp. 
24-25). For infinite-state absorbing chains such as the one we have created 
here, some truncation is necessary to achieve lower and upper bounds. This 
can be done using ideas similar to those used in the proof of Theorem 15.11 
we omit further details. 

The upshot is that is takes an average of 6.07912690331468130722 . . . 
steps to reach coalescence. This is much closer to the lower bound given by 
Theorem 15.11 of 5 than to the upper bound of 15 . Our exact calculations 
confirm results we obtained from simulations. Ten million trials (which took 
only a few minutes to run and tabulate using Mathematica code that was 
effortless to write) gave an estimate of 6.07787 with a standard error of 
0.00184. The algorithm took only a single Markov chain step about 17.4% 
of the time; more than four steps about 47.6% of the time; more than eight 
steps about 23.4% of the time; and more than twenty-seven steps about 1.0% 
of the time. In the ten million trials, the largest number of steps needed 
was 112. 

Similar simulations for small values of /5 led us to conjecture, for some 
constant c near 1, the refinement 

Er(/5) = l + (l + o(l))c/3as/3^0 (5.7) 

of (j5.2p . The expansion (j5.7|) (which, while not surprising, demonstrates ex- 
tremely efficient perfect simulation of Vervaat perpetuities when /3 is small) 
does in fact hold with 



oo 

c = ^ 2-Mn(i 1) 1.016. 

i=l 



(5.8) 
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We will prove next that the right side of ()5.7p provides a lower bound on 
Er(/3). Our proof that it also provides an upper bound uses the same 
absorbing-state approach we used to compute Er(/3) numerically in the 
case /3 = 1 and is rather technical, so we omit it here. 

Using (15.6P and (15. Sp together with the inequality > 1 — x for x > 
we find 



Er(/j)-i ^ 



1 - + 1) 



- 

Application of the dominated convergence theorem to the right side of (|5.9p 
gives the lower bound in (|5.7p . where c = E ln(D_i + 1) is given by the 
series (|5.8p . 

6. Related work 

An unpublished early draft of this paper has been in existence for a num- 
ber of years. Perfect simulation of Vervaat perpetuities has been treated 
independently by Kendall and Thonnes [10]; their approach is quite similar 
(but not identical) to ours. One main difference is their use of the multi-shift 
coupler of Wilson [16] rather than the multigamma coupler. The simulation 
with (in our notation) /3 = 10 reported in their Section 3.5 suggests that, 
unlike ours, their algorithm is reasonably fast even when (3 is large; it would 
be very interesting to have a bound on the expected running time of their 
algorithm to that effect. 

The early draft of this paper considered only the Dickman distribution and 
used an integer-valued dominating chain D with a stationary distribution 
that is Poisson shifted to the right by one unit. Luc Devroye, who attended a 
lecture on the topic by the first author of this paper in 2004, and his student 
Omar Fawzi have very recently improved this approach to such an extent 
that the expected number of Markov chain steps they require to simulate 
from the Dickman distribution is proven to be less than 2.32; see [2]. It is 
not entirely clear that their algorithm is actually faster than ours, despite 
the larger average number 6.08 of steps for ours (recall the penultimate 
paragraph of Section [5]), since it is possible that one step backwards in time 
using their equation (1) takes considerably longer than our simple up/down 
choice (recall step 2 in the algorithm at the end of Section H]) . 

Acknowledgements. We thank the anonymous referee for several helpful 
comments, most notably the suggestion that we try to compute the exact 
value of E T in the case /3 = 1 of the Dickman distribution. Research for the 
second author was carried out in part while affiliated with the Department 
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